{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as npl\n",
    "import scipy.linalg as spl\n",
    "import sklearn.metrics as sklm\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd #for debugging\n",
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "from IPython.core.interactiveshell import InteractiveShell\n",
    "InteractiveShell.ast_node_interactivity = \"all\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.1 Constrained Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Julia cls_solve implementation from https://github.com/VMLS-book/VMLS.jl/blob/master/src/VMLS.jl\n",
    "#function cls_solve(A,b,C,d)\n",
    "#     m, n = size(A)\n",
    "#     p, n = size(C)\n",
    "#     Q, R = qr([A; C])\n",
    "#     Q = Matrix(Q)\n",
    "#     Q1 = Q[1:m,:]\n",
    "#     Q2 = Q[m+1:m+p,:]\n",
    "#     Qtil, Rtil = qr(Q2')\n",
    "#     Qtil = Matrix(Qtil)\n",
    "#     w = Rtil \\ (2*Qtil'*Q1'*b - 2*(Rtil'\\d))\n",
    "#     return R \\ (Q1'*b - Q2'*w/2)\n",
    "# end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 697,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cls_solve(A,b,C,d):\n",
    "    m,n = np.shape(A)\n",
    "    p,n = np.shape(C)\n",
    "    Q,R = npl.qr(np.vstack([A,C]))\n",
    "    Q1 = Q[0:m,:]\n",
    "    Q2 = Q[m:m+p,:]\n",
    "    Qtil, Rtil = npl.qr(Q2.T)\n",
    "    first = np.vstack(np.matmul(np.matmul(2*Qtil.T , Q1.T),b))\n",
    "    denom = first - np.vstack((2*(npl.lstsq(Rtil.T,d)[0])))\n",
    "    w = npl.lstsq(Rtil, denom)[0]\n",
    "    return npl.lstsq(R,np.vstack(np.matmul(Q1.T,b)) - np.matmul(Q2.T , w)/2)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 698,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = 70\n",
    "N = 2*M\n",
    "xleft, xright = np.random.rand(M) - 1, np.random.rand(M)\n",
    "x = np.vstack(np.block([xleft,xright]))\n",
    "y = x**3 - x + .4 /(1+(25*(np.square(x))))+np.vstack(.05*np.random.randn(N))\n",
    "n = 4\n",
    "A = np.vstack([(np.hstack([np.flip(np.vander(xleft,n),1),np.zeros((M,n))])),np.hstack([np.zeros((M,n)), np.flip(np.vander(xright,n),1)])])\n",
    "b = y\n",
    "one = np.block([1,np.zeros((1,n-1)), -1,np.zeros((1,n-1))])\n",
    "two = np.block([0,1,np.zeros((1,n-2)), 0, -1, np.zeros((1,n-2))])\n",
    "C = np.vstack([one,two])\n",
    "d = np.zeros((2,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 699,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    }
   ],
   "source": [
    "theta = cls_solve(A,b,C,d)\n",
    "Npl = 200\n",
    "xpl_left = np.linspace(-1,0,Npl) #continuity constraint 1\n",
    "ypl_left = np.matmul(np.flip(np.vander(xpl_left,4),1) , theta[0:n])\n",
    "xpl_right = np.linspace(0,1,Npl) #continuity constraint 2\n",
    "ypl_right = np.matmul(np.flip(np.vander(xpl_right,4),1),theta[n:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 700,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x12704bda0>"
      ]
     },
     "execution_count": 700,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12700ab00>]"
      ]
     },
     "execution_count": 700,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1270304a8>]"
      ]
     },
     "execution_count": 700,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAFpCAYAAAB0yyjhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNXd+PHPmclkYwtL2MKOEERZAgFUFAUUcMeliFWrVX/WWp+u8hSrVmsXqdTHPrY+WqvWta4oooKIAlYRUCDssgtCwg5hS0ImM+f3x8yEyeTOlrkzd5bv+/XileTOnbnf3Azfe+bcc75Haa0RQgiRWWxWByCEECLxJPkLIUQGkuQvhBAZSJK/EEJkIEn+QgiRgST5CyFEBpLkL4QQGUiSvxBCZCBJ/kIIkYEk+QshRAbKsjqAYNq1a6d79OhhdRhCCJFSli9ffkBrXRhuv6RN/j169GDZsmVWhyGEEClFKbUjkv2k20cIITKQJH8hhMhAkvyFECIDSfIXQogMJMlfCCEykCR/IYTIQJL8hRAiA0nyF0KIDCTJXwghMpAkfyGEyECS/IUQIgNJ8hdCiAyUtIXdhEgWM8vKmT53IxWV1XQuyGPK+GImlhRZHZYQMZHkL0QIM8vKufedNVQ7XQCUV1Zz7ztrAOQCIFKadPsIEcL0uRvrE79PtdPF9LkbLYpICHNIyz/FSZdEfFVUVke1XYhUIck/hUmXhDEzL4idC/IoN0j0nQvyYg1TCEtJt08Ki1eXxMyyckZOm0/PqR8yctp8ZpaVx/R6ieS7IJZXVqM5dUFs6u8wZXwxeQ57g215DjtTxhebEK0Q1pGWfwqLR5dEsE8Ty3YcYsGG/RG3pq3qjgp1QWzK8X3Pka41kW4k+aewYF0SrfIcjJw2v0nJKljyfHXJd2jvz+G6l6zsjorHBXFiSZEke5F2pNsnhRl1SThsihO1dU3u9giWJHXAz6G6l6wcIROsLz6aPvpU7vYSIlKS/FPYxJIiHrl6AEUFeSigqCCP5rlZOF0NU3U0iTeaJBltKzsRI2Ri7aM3+56BEMlKkn+Km1hSxKKpY/h22qUsmjqGyiqn4X6RJl6j5KmC7BttKzsRI2SMLoiPXD0g5m4vGdcv0o30+aeZWIcmGt3gHN2vkBnLyxskxVCt6Snjixv0+Yfb32zR9NH735guyHdwOMjFs7yyusn3UYRIRpL804wZidcoeZZ2bxPxiJdUGSETeGM6WOIHz6cf30VV5lOIdKC0DryVlxxKS0v1smXLrA4jJUUyzNKsoZjHT9axee8xvjtUxcHjtRw8cZITJz3J1KYUOQ4bhc1z6NAyl04FufTr2IL87ORoc4ycNt/wU1IgReMb3uDpUlo0dYzpcQkRC6XUcq11abj9kuN/oTBVuG6Ppg7F1Fqz42AVX249yKKtB1i1s5JdhxsmT5uC5jlZnmSpoabO1eAGtFLQu7A5A4taMapvIef3LaR1s+ym/7IxiPQ+SLDmkZR4EKlMkn8GinYi1JZ9x5m1qoIPVlWw7cAJADq0zKG0RxsmD+tK3w4t6NmuGe2a59Aqz4HNduoWsdaaw1VO9h2r4buDVazffZS15Uf5bNN+3ikrx6agpFtrrhnShcsHdaJFriO+v7yfYPdHAtmVwmXwCVlKPIhUJsk/A0UyFPNknYvZa3bzwqLtrNp1BKXg7F5tuWVkD0ae1o5e7ZqhVLBxQKcopWjTLJs2zbLp17El487oCIDbrVlTfoT5G/YxZ+1ufvPuGn7/wXouHdiJO0b1om+HFub8siEY3R8x4tKaPIfdshvYQsSDJP8k498X3yrPgVJQWeU09aZpqBFBx2qc/GvRdl5avJ0Dx2vpXdiMBy7rz+UDO9G+ZW7Mx/ax2RSDuhYwqGsBP7+wD6t2HeGNr3fy3spyZqzYxcVnduTu0X3o37llzMcKdn8j8Ma0UuA26OMp8j4n2W9gCxENueGbRAL74gPlOexRjVmP5ji5WTbGnt6BRVsPUFnlZHRxIbee25NzT2sXUQvfLIdO1PL8F9/y4pfbOXayjmuGdOHXFxfTvkXTLjxGv2uw8xjNvkIkq0hv+EryTyKRjD4xa4SJrzVcXllN63wHLrfmaE0do4sL+eVFxQzo0irmYwQeK5pW85FqJ08t3MpzX2wjJ8vOzy/swy3n9CDLHt28xGDnNNh5lPURRKqT0T4pKJLRI2aNMJlYUsRp7Zvz4Kx1LN9xmEFdC/jtZf0Z2r21Ka/v09SRRa3yHEy9uB+TSrvw8Afr+cOH3/DB6t08ft1gerZrFvHxoy01IUXcRKaQ8g5JJJLRI2aMMDlZ5+Ivczdy5ZOL2HHwBI9eO5B3f3yO6Ykfgo8s+t376yIqntarsDn/umUYT1xfwrb9x7nkfz/nlSU7iPQTq5WlJoRIZtLyTyLhRp+EG2ESrMvCf3u75jlk2RW7j9Rw7dAuPHBZf1rlxW94ZbAW9uEqZ/2M2nCfBpRSXDGoM8N7tGHK26u4f+ZaFm87yPRrB4adMBbJjGfp6hGZSPr8k0xTRvvMLCvnoVnrqKxuWJ4gz2HnmqFFjeryANxxXi9+c+npcf1dIPJZtBDZ/Qy3W/PM59t49KMN9O3QgmuHduFfi7Y3eTaz3OQV6UZu+GaIcCOEgk1QSlRpgnDx+VPAt9Mujeh1/7NpP3e+spyq2oavG03inllWzq/eXGXp+RHCbJEmf+nzT3FGfer+jBIbJK40gVGJ5YIg3UzR9MOP6ltIi5zGXT6Rll/2XZSsPj9CWMWUPn+l1ATgfwE78KzWelqQ/a4F3gKGaa2lWW8g0v5n/6GaTZHIG56BI2iCdbVEO2N237GThtsjSdzhLpoF+YkrMyGEFWJO/kopO/AkcBGwC/haKTVLa70+YL8WwE+BpbEeM11FMixyZlk5v3t/Xcjyw/7aNc/meE0dNXXu+m1NLU0Q7sIU6YXLrJLPsaxdEO4Ccbymjpll5dLvL9KWGd0+w4EtWuttWuta4HXgSoP9fg88CtSYcMy0FG4VKd/FIdLEX9q9NV/8egzTrhnY5JWtfMItbxjt8oeBK5A1JckarToGnt87nHAXCKdby+pdIq2Z0e1TBOz0+3kXMMJ/B6VUCdBVa/2BUuoeE46ZlsJNSArXVQGQk2XjZJ2bKeOLueuC3iilTJm4FK4SaKjx/PFqPQd+gujUKpf2LXN5b1UFg7oWcOu5PYM+N5KibtLvL9KZGcnfqPBL/V00pZQNeBy4JewLKXUHcAdAt27dTAgttYTrxgiXjBx2hdPl5tFrBjJpWFdTYwt3YQo1nt+o+8SssfWBFzany81//buMhz9YT47Dxg0jugd9HhDyvolMBBPpzIxun12Af6bpAlT4/dwCOBNYqJTaDpwFzFJKNRqKpLV+RmtdqrUuLSwsNCG01GLUjeHfPx8qGSlAa3jqxqGmJ/5Qx/ZtDxVbYPdJtF1E0XDYbTxxfQmjiwt5YOZa/vDB+qAziX1dT3+9bnDI8y5EOjIj+X8N9FFK9VRKZQOTgVm+B7XWR7TW7bTWPbTWPYAlwBUy2qcxo2GR/v3zwfq4bQqys2y8fNsIxnvr5TfFzLLyoIky3IUpVKIM/FQQ7t5GrLKzbDx5wxC6tM7n2S++DXuRCXfeA4U6T0Kkipi7fbTWdUqpu4G5eIZ6Pq+1XqeUehhYprWeFfoVhL9Q/fOBfdwdWuZS53Zz0unm5dtHMLhrQZOPG26kUbgROhNLigxnGUPjTwXRFltrivzsLGr9Rjj5BFuxLNL7Ik0tVCdEspEZvilqz5EarntmMYeO1/LSbcMp6RZbUbZoSx8bCVUqAU5dOGwJmnXcc+qHhuvvRjOTOFDJwx8bjraSGcEiWUhJ5zS292gNk59ZzEGTEj+Y0xoP9ukAaHBRMEr88ehjj2UegJGZZeVBh9nKyCCRaiT5p5jKqlpuem4p+4+d5KXbRjDEhMQP5iVKo+6TkdPmGw6ptCuFW+u4VdI0Gs6ZZVNNvsiEuichI4NEqpHkn0Kqauu49YWv2X6gihduHWZq/f1ISh+H4j90syDfgdZQWe0MWlgOwK11k7pfAmc5F+Q5eOiKMwz78eHUcM78bDtVtS7yshvfNI9EqNa9jAwSqUYKu6UIp8vNXa+uYOXOSp64fjDn9G5n6utHO+LFX+DQzcNVzvobv8ESPzSttTyzrJwpb69q0P1SWe1kylurDEfdTCwpYsr4YooK8qiqdeGwK372ehk7D1VFfexg8RbkOeRmr0g50vJPAW63Zspbq1i4cT+PXD2ACWd2istxmjITOFRZ5FCa2sc/fe5GnK7Gx/KVYwi3KLvTpXG6NDc+t5RPfnk+Du+awJFMOgv26eihK86I+vcQwmrS8k8B0z/eyMyVFUwZX8z1w5Nn5nO4ssjBNLW+EITuejF6LFhJjB0Hq/jLxw1rJpk9H0CIZCYt/yT35tc7eWrhVr4/oht3XdDb6nAaiKTWUKBYh0QGuzHteyxQqIvFPz7bxtm92oatW+RPFngX6UJa/kls0ZYD/ObdNZzXpx2/u+IMlDIqo2SdaIc3mjGcc8r4Yhz2xufBEWQUT/B++iyKO7Rgyturg15MZPimSGeS/GMUr6n+W/Yd485XltOrsBlP3jCkvm86mQRLrDZF/Wpddu8FK9IuknDnc2JJEdOvHURrv8VWCvIcTP/eIMPXDnaxOHHSxeWDOnH4RK1hyYxQv58Q6UBm+MYg0hmt0Y5jP3yiliufXERVrYuZPzmHLq3zTYvXjEqa/q9n5uLn8VpMffDvPjYsO1FUkMd1w7ryP/M2kW23UetquOCN0XHNPodCmE1m+CZAqBr2NU53k+q/uNyan75exp4jNbzxo7NMTfxm16Qxa0Uun2j63qNxxCDxg6db58cX9Gbe+r1sO3CcNs1y2Xu0JujvIXV9RDqR5B+DUDXsA0WaxKbP3cjnmw/w52sGmFK2wf9145FYmzo81OiCEa+Cb6FmLzvsNh6bNIjL/vYFg7sW8PRNQ4O+TrBz+Ks3VwFyARCpRZJ/DEKNPDESLol9uHo3T3+2lRtGdOO6YeYO6UxEJc1IhGo9m12Lxyfc7OW+HVrws7F9mD53I/PW7+Wi/h0MXyfYuXJpzb3vrGHZjkMs2LBfuoRESki+u4gpJFh9/WBCJbGNe44x5e1VDOlWwIOXmz9pKNxiLIkS6hNIuDUDmiqS8fl3jOpFcYcW/Pa9tRw/WWf4OqHOVbXTxatLvovLAjVCxIMk/xj4J5VwQiWxI1VO7nh5Gc1ysnjqxqFkZ5n/Z4lXYo1WqE8g8ZxEFW7BeIfdxp+uHsCeozU89rFxAbdwF/vAoRNmLlAjhNmk2ydGvj7vYPXwwZPEgnUBaK351VurKD9czet3nEWHlrlxixPMuzkbDf8+/mC1/H2taisnUQ3t3pobRnTjxS+3c1VJEQO7NFwcxxdXNOUsZK6ASFaS/E0SrF85XMv1+UXb+eSbvTxwWX9Ke7SJa4xWJNbAPv5E1fJvqv+e0I+P1+3lvnfX8t5PRmKzNZwj4Dt/gX9rReOWP8hcAZG8pNvHJE3psli1s5Jpc77hov4duHVkj4TFmkjBSkDYlUrK+jgtcx3cd+nprCk/wlvLdxruY/S3vuGsbknRrSZEpKTlb5JoJ/8crXFy92sraN8il+nXDky60g1mCdbt0dRa/olwxaDOvLx4B49+tJEJZ3aiVZ6j0T5Gn6JKu7eRCWAiZUjyN0G0k3+01kydsZqKyhre/NHZFORnJzTeRIrX8M14Ukrx0BVncPnfv+CJTzfzwGX9I3qeFH0TqUS6fUwQaviikVeXfsfsNXuYMr7Y1NW4klGyjDKK1plFrZg8zHPzd8u+Y1aHI4TpJPmbIJoJVBv2HOXhD9Zzft9C7jivV7xDs1wq18C/Z1xf8rLtPPzBN1aHIoTppNvHBJF2bZysc/Hz11fSMjeLxyYNajSSJF2landI2+Y5/GxsH/7w4Td8vnk/5/UptDokIUwjLX8TRNq18T8fb2LDnmM8eu1A2jXPSWSIooluOrs7XVrnMW3OBtzu5KyAK0RTSPI3QbiujZll5ZT+YR7/+M828rPtHK02Lh8gkk9Olucivq7iKLNWVVgdjhCmkXr+cTazrJypM1ZTUxe+VrxITm635oonv+DwCSef/up8cqOo5yREokVaz19a/nE2fe7GBokfpOZLU8RrxbRI2GyKey8+nfLKal5ZsiNhxxUinuSGb5zJ+rCxS4ZFVEae1o7z+xby9wVbuG5YV1rkNp74FY6sAiaSibT842j/sZMEG9CTzJOckk208yji5Z5xxVRWOXn+i+1RP9d3AZOSzyJZSPKPowdnrUUpRU5AieZUmOSUTJJlIZoBXVox/owOPPv5NiqraqN6brJcwITwkeQfJ7PX7Gb2mj388qK+/PmagSk5ySlZJMtCNAC/uKgvx2vreH7B+qielywXMCF8pM8/Dg6dqOWBmWsZUNSKH43qRZbdJsk+BuGWYUykfh1b8nD3NVz+1f/j4LCvads+sr9rKtY4EulNWv5x8NCsdRytcTL9ewPJssspjoXvJmm104XdW/nU6k9PF4yZQEtVRd6ypyN+TqrWOBLpSzKTyT5et4dZqyq4e3Qf+nVsaXU4Kc3/Jil4FoLxJUwrP0l17TsY2xkTyV/5PFQdiug5qVzjSKQn6fYx0ZEqJ/fNXMvpnVpy1+jeVoeT8kLdJLU8aZ53D6x7F756Bi6YGtFTUrXGkUhP0vI30SNzvuHQiVqmXzsQh3T3xCypb5J2PBOKL4Ul/wc1R62ORoioSYYyyVffHuL1r3dy+7k9ObOoldXhpIVkGuVjaNQ9UHMEvn7W6kiEiJokfxPU1rn5zbtrKCrI42cX9rE6nLSR9DdJi4bAaRfC4r9D7QmroxEiKpL8Y+Sr2Lll33FqnC4+XrfX6pDSRkrcJB3131B1EJa/YHUkQkRFqnrGYGZZOb+esZqTUrEzs71wGRzYDD9fDVmyToOwllT1jEGkFSQf/WhDg8QPMmU/I533Szi+B9a8ZXUkQkRMhnp6+SYTlVdWowDf56FQFSQrjtQYvlZSjEYRidNrNHQcAIuegEHfB5u0qUTyk3cpjScTBXaEGbXmD5+olYqdwkMpOOencGAjbP7YcBcr1yMQwogkf4wnEwUKbM3/+aMNAFKxU3iccRW07AJfPtHoISnnLJKRJH8i66bxb82v3FnJG8t2ctu5PaVip/CwO+Dsu2DHIti1vMFDUs5ZJCPp8yd4xUUf/9a826158L21FDbP4adj+9Ai1yHJXngM+QEs/DN8+b8w6aX6zUk9U1lkLGn5YzyZyNed79+an1lWzpDfz2PVriMcr6lj1KMLpA9XnJLTAobdBt+8Dwe31m9O+pnKIiNJ8sd4MtHj1w1m+7RLWTR1TH3inzpjNZXVTgCqnC4OVzmlD1c0NOJHYMuCxU/Wb0r6mcoiI0m3j1e4iovT526kJmBMv7+kqTYprNWiIwycBKteg7EPQF7r+veELN4ukokk/wiFuifgI324AoARP4ayV2DFSzDyZ4CUcxbJx5RuH6XUBKXURqXUFqVUo+LmSqlfKqXWK6VWK6U+VUp1N+O4iaK1JjuCEs3ShysAT7nnHufBV/8EV53V0QhhKObkr5SyA08CFwP9geuVUv0DdisDSrXWA4G3gUdjPW4izVxZTq3LjcMeZFYX0ocrAoy4E47shI0fWh2JEIbMaPkPB7ZorbdprWuB14Er/XfQWi/QWld5f1wCdDHhuDGLZNbl8ZN1/Gn2BgZ1LWgwpr8gz0HrfIeM7xfGii+Ggu6wJPJ1foVIJDP6/IuAnX4/7wJGhNj/NmCOCceNiW/WpW/yTbAaPk8t3ML+Yyd55qahlHRrzdVDkuK6JZKdzQ7D74CP74Pdq6DTIKsjEqIBM1r+Rn0hhnWilVI3AqXA9CCP36GUWqaUWrZ//34TQgsuklmXuw5X8c/Pv2Xi4M6UdGsd13hEGiq5ERzNpPUvkpIZyX8X0NXv5y5AReBOSqkLgfuAK7TWJ41eSGv9jNa6VGtdWlhYaEJowUUy6/LRjzaigCkT+sU1FpGm8gpg8Pdh7dtwfJ/V0QjRgBnJ/2ugj1Kqp1IqG5gMzPLfQSlVAvwDT+JPiv8F4WZdrvjuMLNWVXDHqF4UySge0VQjfgSuWlj2L6sjEaKBmJO/1roOuBuYC3wDvKm1XqeUelgpdYV3t+lAc+AtpdRKpdSsIC+XMKFmXWqt+f0H62nfIoc7z+9tUYQiLbTrA33GeRZ5rzP8wCuEJUyZ5KW1ng3MDtj2W7/vLzTjOGYKNevyvZXllH1XyaPXDqRZjsyDEzEacSe8cjWsmwmDrrM6GiGADJ/hazTrssbp4s9zNnBG55ZcKyN7hBl6jYa2p8Gy5yT5i6Qhhd0CPPv5NiqO1HD/pf2xBVuqS4ho2GxQeivsXAp71lgdjRCAJP8G9h2r4f8WbmVc/w6c3but1eGIdDL4+5CVB18/Z3UkQgCS/Bt4bO4mnC43v7nkdKtDEekmrzUMuAZWvwk1R6yORghJ/j7f7D7Km8t3cvPZPejRrpnV4Yh0NOx2cJ6AVW9YHYkQkvx9Hv1oAy1ysrh7zGlWhyLSVecSKBrqGfapDSfBC5EwkvyBxVsPsmDjfn4y+jQK8rOtDkeks9Lb4MBG2P6F1ZGIDJfxyV9rzbSPNtCpVS43n9PD6nBEujvzasgt8LT+hbBQxif/OWv3sGpnJb+4qC+5ATN+hTCdI89T8G3DB3Bsj9XRiAyW0cnf6XIzfe5G+nZozjUyoUskSumt4K7zLPMohEUyOvm/8fVOvj1wgl9P6IddJnSJRGnbG3qP9RR7k2UehUUyNvmfOFnHXz/ZzPAebRjTr73V4YhMM+w2OFYBmyxf10hkqLSu7TOzrNywcBvA8198y4HjJ/nHTUNRSlr9IsH6jIeWXeDr55hZMyTo+1SIeEnb5B9qmcbz+rTjH//ZxoQzOjK0u6zQJSxgz4IhP4CFf+Lvmz+l3OkpJxJsOVEhzJa23T6hlmn82/wtVDtdTJlQXP9YJIu5C2GqkhtxYWOi/qTB5sDlRIWIh7RN/sGWaSyvrObVpTuYVNqV3oXNgVOfEsorq9Gcan3JBUDEVasiFroG8T37Z2TR8MZvsPevEGZJ2+QfbJnGfIcdm1L8/MI+9dsiWcxdiHiYm3sxHVQlY2xlDbYHe/8KYZa0Tf5GyzTmZNmornPxg7O706Flbv32SBZzFyIeRk6YzF7dmuvt8+u3+ZYTFSKe0jb5Tywp4pGrB1BUkIcCigryOL1TS/Id9kbr8oZbzF2IeLlyaHcOF0/mfPtqijhAUUEej1w9QG72irhL29E+0HCZxnUVR7j0iS/4rzGn0bZ5ToP9powvbjAyCKT1JRKn3yV3waanWTS+HEbfbHU4IkOkbcs/0OPzNtEyN4vbz+vV6DGjTwnS+hIJU9ANThsLK16WGb8iYdK65e+z4rvDfPLNPqaML6ZVnsNwH6PF3IVImKG3wBs3wpZPoHiC1dGIDJARLf//+XgTbZtlc4uUbBbJqu8EaN4Blr9gdSQiQ6R98l+89SBfbDnAjy/oTbOcjPigI1KR3QGDb4DNc+FohdXRiAyQ1slfa81jH2+kY8tcbjyru9XhCBHakB+AdkPZK1ZHIjJAWif/hZv2s2zHYe4ec5os1CKSX5ue0Gu0p86/2xV+f5GeZv83zLwr7odJ2+Tva/V3bZPHpNKuVocjRGSG3gxHdsLWBVZHIqzy7X+g6mDcD5O2yX/uuj2sLT/Kz8b2Zfaa3VK0TaSG4kshvx0s/5fVkQgr1J2Eg5uhwxlxP1RaJn+XW/PYx5voXdgMG0jRNpE6srKh5AbYOEfW+M1EBzZ5lvhs3z/uh0rL5L/vWA1KwS8u6stj8zZJ0TaRWobcDNoFK1+1OhKRaHvXe752ODPuh0rL5N+pVR5zfjaKS87sFLK0sxBJqW1v6H6uZ9SP1lZHIxJp71qwZ3veA3GWlskfwG5T2GwqaHE2BdL1I5LXkJvg0DbYscjqSEQi7VsPhcWeeR9xlrbJ32fK+GKMVujVIF0/InmdfgXktPTU+xGZY+86aB//m72QAcl/YkkRwT44S71+kbSy82HAtbD+Pag5YnU0IhGqDsGx3QkZ6QMZkPzBU6XTiNTrF0mt5Caoq4Y1b1sdiUiEves8XzvEf6QPpGlVz5ll5Uyfu5GKymo6F+Qxul8hM5aXS71+kVo6l3hGfZS9DMNuszoaEW/7EjfSB9Kw5W+0GPuM5eVcM7RI6vWL1KKUp/VfUQZ71lodjYi3vWshr42numsCpF3LP9hi7As27GfR1DEWRSVEEw2cBPMe8LT+L/6z1dGIeNq73tPfr4yGqJgv7Vr+shi7SCv5baDfZbD6Dc/Uf5Ge3G7Y903CbvZCGiZ/WYxdpJ0hN0H1YdjwgdWRiHg5tA2cJxLW3w9pmPynjC8mL6B8s9zcFSmt5wXQqpuM+U9ne1Z5vnYamLBDpl3yl8XYRdqx2TzF3rYthMrvrI5GxMOeNWBzQOHpCTtk2t3wBVmMXaShwTfAwmlQ9iqMvtfqaITZdq+Gwn6eqq4JknYtfyHSUkFX6D3aU+lTVvlKL1rDntUJ7fIBSf5CpI6SmzyrfG1baHUkwkzH98KJ/dAxsck/Lbt9/AXO9p0yvli6hERq6nepZxLQipfgtLFWRyPMsnu152vHAQk9bFq3/I1m+8oqXiJlZeXAwOtgw4dwIv5rvIoE8Y306Zi4YZ6Q5sk/2GxfKeUsUtaQm8Dt9Ez6Eulhzxpo3QNyWyX0sGmd/GW2r0g7Hc6AzkM85R5kla/0sHt1wvv7waTkr5SaoJTaqJTaopSaavB4jlLqDe/jS5VSPcw4bjgy21ekpSE3eSpAlq+wOhIRq5ojcPjbhI/0AROSv1LKDjwJXAz0B65XSgUWpL4NOKy1Pg14HEhIhSqZ7SvS0pnXQFYelL1kdSQiVr4a/ina8h8ObNFab9Na1wKvA1cG7HMl8KL3+7eBsUrFv3SdzPYV6WRmWTkjp82n50NfMMc9Aueqt6D2hNVhiVjUj/RJfPI3Y6hnEbDT7+ddwIhg+2it65RSR4C2wAETjh+SzPYV6cA3cs03gOEaZdXIAAAgAElEQVRf1edxcc5Cls95gaFX/sTi6EST7VkDzQqhRceEH9qMln+w9dGj3Qel1B1KqWVKqWX79+83ITQh0kPgyLWvdD+2uTtiX/mKhVGJWFV+u4yl1UX0vHc2I6fNT+gwdDOS/y6gq9/PXYCKYPsopbKAVsChwBfSWj+jtS7VWpcWFhaaEJoQ6aHxCDXFW64LGKzXw4EtlsQkYvPe8m/Jr9zMitqulsxDMiP5fw30UUr1VEplA5OBWQH7zAJu9n5/LTBfaxmnJkSkjEaove06jzptY9NHT1kQkYjVrLnzyFYu1rp71m9L5DykmJO/1roOuBuYC3wDvKm1XqeUelgpdYV3t+eAtkqpLcAvgUbDQYUQwU0ZX9yo73Q/rVngHkzrLW+Dq86SuETTdTzxDQCrdK8G2xM1D8mUcf5a69la675a695a6z96t/1Waz3L+32N1vp7WuvTtNbDtdbbzDiuEJliYklR45tkwJuuCyikkiVzX0t4TCI2Z+Xu4JBuzi7dsIs7UfOQ0nqGrxDppMggKSxwD2a/boVz2YsGzxDJbFSznazlNPzHwyRyHpIkfyFShFFSqCOLGa5RnO1aDsf2WBCVaJLaE7Q6toX2xWdbNg8p7Us6C5EuJpYU8dCsdVRWOxtsf9N1PndmvQ+rXoNzf2FRdCIqu1eDdtNvyCgWfX+MJSGkbfKXOv4iHQS+jy8b1IkZy8sbjPnfndWVA22G0K7sFRj5c4j/5HkRq4oyz9eiIZaFkJbdPlLHX6QDo/fxjOXlXDO0qFFXQbvzboeDW+C7xVaHLSJRsQJadLZkZq9PWrb8Q9Xxl9a/SBXB3scLNuxn0dSAroLaiTDn17DiZeh+TgKjFE1SvsLSVj+kactf6viLdBDV+zi7GZx5NayfCTVH4xyZiEl1JRzaCp0HA34F+6Z+mNASD2mZ/KWOv0gHUb+Ph/wAnFWw9u04RiVitnul52vnIZZ2Uadl8pc6/iIdRP0+LhoK7c+A5TLmP6n5FuHpXGLpUrNp2efv69eX0T4ilUX9PlbK0/r/6NewexV0GtTgYRkBlyQqVkDrnpDfxtIu6rRM/iB1/EV6iPR97EvsJypbszTXQfncp+h1y9MNHvdfD8DXveA7hkig8jLoOhzwdOGVGyT6RHRRp2W3jxCZxL/fuJLmzHYNp9237/H+16dKPVvZvSD8HN8HR3fVj/Sxsotakr8QKWxmWTm/enNVg8T+et0YWqoqFr77z/rRIzICLkns/MrztcswwNqlZtO220eIdOdr8bsClsZYqvux1d2JyVnzmVE5invfWUNBvoPDVc5GryEj4BJs11dgc0CnwfWbrOqilpa/ECnKqCvHQ/GG6wKG2TZxmtpFtdOF1sgIuGSw8yvPjXhHrtWRSPIXIlWF6rKZ4RpFrbYz2b4AgCPVTsu6F4RXXa2npk/XEVZHAki3jxApxX+4pk2pRl0+PgdpxTz3UK62f86jdZMpLGgpI+CstncN1NVA12FWRwJIy1+IlBE4GzRY4vd53TWGNuo4l2cvl+6dZFB/s3e4tXF4SfIXIkUE7+NvzK4Ui9xnUqHaM7X9UmnxJ4OdX0HLLtAqOf4W0u0jRIqIZlimW2u+nXY5fLYBFvwBDm2DNr3CP1HELOhM6p1f1U/uSgbS8hciRUQzLLN+35IbQNlgxUtxikr4C1ao7aMvV3gmd0nyF0JEy2g2qMOucNgartzVYAhny87QZzys/De4Go/zF+YKNpP68wWzPT9I8hdCRMtoNuj0awcx/XuDQg/hHHozHN8Lm+ZaFXrGMKrTA9Czei1k5UKHAQmOKDjp8xcihQQbrunb5utv/sUbK0/1Nw+8CFp0ghUvwumXJTrkjGIPMvx2iG0LdB4CWdkWRGVMWv5CpIlg/c33v/8NL1Sdi2vTPK565E1ZyzqOjBJ/DrWcqbYlzfh+H0n+QqSJYP3Nry75jmerzkUBo45/lLCVojJRkcFN+YFqG9nKlTQze30k+QuRJoINBdXALl3IF+4zmZS1kFpnrZRyjhOjm/IjHRvRKOh2tkVRGZPkL0SaCDcU9FXXWIrUQUbbVkop5zgxuin//Q47UR3OgPw2VofXgNzwFSJNTBlf3GC1LgCFp+UP8Il7KHt0a260f8I3Lc+1JMZM0OCmvMsJ026EkhutDcqAtPyFSBOBrc7W+Q5yHaf+i7uw87prNKNsq3nw3HzrAs0kFSvBeQK6n2N1JI1I8hcijUwsKWLR1DE8ft1gapxuqp3uBo/PdowDm41x1bMtijDD7PjC87X7SGvjMCDJX4g0FKwI3ImcDtiKL4ayV6DupAWRZZjti6BdMTQvtDqSRiT5C5GGQq7ZO+w2qDoI699LcFQZxlUH3y2BHsnX6gdJ/kKkpWAjfzoX5EHPCzwVPr9+LrFBZZo9q6H2WFJ2+YAkfyHS0uh+haiAbfUF32w2KL0Vdi6BvessiS8j7Fjk+dojOUdWSfIXIs3MLCtnxvJy/AsNKOCaoX5DEAffAPYcaf3H0/ZF0KY3tOhodSSGJPkLkWaMbvZq4LWlO0+VdchvA2deDavfgJPHEh9kunO74Lsvk7a/HyT5C5F2gt3sdWndsK5P6W1QexxWv5nA6DLE3rVQcyRp+/tBkr8QaSdUmYdqp+tUXZ8updBxACx7HsIsBi+itO0zz9eeo6yNIwRJ/kKkGaPiYv4qKquZWVbOyD8v4N6dI2DvWv7z6fsJjDADbFsAhf08K6klKUn+QqQZX5kHuwoc7+NRkO+or/s/03UOR3Q+J/7zf1Lm2SzOGtixGHpdYHUkIUnyFyINTSwp4rFJgxp9Ashz2NGa+hvC1eTyums0F6mlvPDRIitCTT87l0JdNfQaXb9pZlk5I6fNp+fUDxk5bX5SXGgl+QuRpozKCz9y9QCOVDdcyP1l1zgUmgtPfGhNoOlm20JQ9vqRPsFWWLP6AiAlnYVIY0Zr/k6fu7HBQuO7dCGfuIdyQ9Z8T5eFIzfRYaYU3zrJFZXVp9ZJ9j/H2xZCl2GQ0wIIvsLa9LkbDddjThRp+QuRYYxuCL/GxbTmKKydYVFUqSFsK77qEFSUQe9TXT4h6yxZSJK/EBnGqDto4lWTofB0WPq0DPsMIVQrHoDtnwO6wc3ekHWWLCTdPkJkIKPuINw/gg9+7qlE2T251ptNFmFb8Vs+hewWUDS0/jGjFdbq6yxZSJK/EGkobL+0wf5/m9eBGboZy1/8Pccu/6el/dHJqnNBXoP7Jf7b0Ro2z4Peo5m5el+D83/N0CIWbNgf8d8jEST5C5FmfP3Svpamr18aMEw4p/bXvJ41mtuZzYXvLAQusDxBJZuQrfi96+BYBStyhjU6/zOWl/PI1QOS6nxKn78QaSZsv3SI/V9xXYRC8z39UdD9M1mw4bMTS4pg88cA/O6bzlGdf6vE1PJXSrUB3gB6ANuBSVrrwwH7DAaeAloCLuCPWus3YjmuECK4UP3SRt1BFQHDPue5S/m+fT5PVk5MVMgpxfB+CXi6fDoOZPV24xu5Vo/uCRRry38q8KnWug/wqffnQFXAD7TWZwATgL8qpQpiPK4QIohgo0ha5TkMhym2ynM02O+fdZfQWh3ntuaLExBtmqiu9Mzs7TMuaUf3BIo1+V8JvOj9/kWgUVNBa71Ja73Z+30FsA9IvtWMhUgTRuP48xx2lMKwO0IpGuy/XBdTpvtyZ/ZHnrr0IrxtC0C7oM+4oOff6tE9gWJN/h201rsBvF/bh9pZKTUcyAa2xnhcIUQQwfqlK6uchvtXVjkb7V87/Cc0q9oJGz5IaOwpa/M8yGsNXUpD3xdIIkqHmdChlPoEMFqH7D7gRa11gd++h7XWrYO8TidgIXCz1npJkH3uAO4A6Nat29AdO3ZE8jsIISIwctp8w2GKRQV5LJo6puFGtwv+NhTy28Ltn0CQCqECz7n6Sx/oPQauedbqaFBKLddal4bbL2zLX2t9odb6TIN/7wF7vUndl9z3BQmmJfAhcH+wxO891jNa61KtdWlhofQMCWGmqLojbHY4+ydQvswz6UsEt3MpVB2EfpdaHUlUYu32mQXc7P3+ZuC9wB2UUtnAu8BLWuu3YjyeEKKJfN0RrfNP3eDNyQqRAgbfAHlt4Mu/JSC6FLbhQ7Bnw2kXWh1JVGJN/tOAi5RSm4GLvD+jlCpVSvk+/0wCRgG3KKVWev8NjvG4QogmqnG667+vrHby8zdWUvLwx41LDGfnw7DbYeNsOLA5wVGmCK0990V6XVBfxTNVxDTOX2t9EBhrsH0ZcLv3+1eAV2I5jhDCHEYTwAAOVzmNZwEP/3+w6H9h8ZNw+V8TFWbq2LsODm+Hc39Rvyna0hpWkRm+QmSQUBONDGehNm8PgybDyn/D8f1xji4FbfgQUFB8CZC8C7cYkeQvRAYJN9HI8OJwzn+BqxaWPpWUyxFaasMH0HWE5yJJ9KU1rCTJX4gMYjTix5/hxaFdH+h/Jc7F/+BP7yxJiVZtU0V1cTu8Hfashn6X1G9K1oVbjEhVTyEyiK/v+aFZ66gMWMs35CzUUffgWD+TSe45/J2r6jcnw3KETRXYNz+6XyEzlpdHXA2Vte94vvY/VdggZMnnJCMtfyEyzMSSIlY+OI6/Xjc48lmoHQfwqauE27LmkE9Ng4eSsVUbjlHf/KtLvouuy2bdO561elt3r9+UKqUdQFr+QmSsoNUpg3gj7zrG1k7lBvsn/NN1Wf32ZGzVhmPUNx+s1oH/xc33aSH3yFY+zVnDmjOnMsBvX9/5TIXRPpL8hRARueTiK/hy5qvckfUhL7nGcZLspG3VhhPNpxXfxc1/kZyf2pfg1oq7V3XnF73LGyT3aC+qVpFuHyFERCaWFOE69x4K1REm2xckbcGySAT7tBJYwcj/4nbq04LmcvtivtbF7HC2SsqRPJGQ5C+EiNh5F14J3c7md20/ZdE956Zk4ofgffM3nNUt6H0Q36eFYrWTPrZy3ned3WB7qpFuHyFE5JSCUVPglathxYueGcApqCl9876RPFfZF1GnbcxxDa/fnook+QshotN7DHQ7B/7zF0/xt+x8qyNqkmj75qeML+b+d1Zyle1zFrhLOEirlL3nAdLtI4SIllIw9gE4vge+tr5+faJMLCnimZFH6aAqmeE6L6XveYC0/IUQTdH9HOg9Fr54HIbeArktrY4opKYWWwt83htt34O8Njx9/32QlZ2AyONHkr8QomnG3Af/HANLnoILfm11NEH5D9GE0DN3/ZN9Qb6D4zV1ON2eGQDHKvdTWP0JW3tNoneKJ36Qbh8hRFMVDYV+l8Hiv0PVIaujCSpYsbVfvbmqQe2ewFm/h6uc9Ykf4HL7EnKUkz+UlyQq9LiS5C+EaCTiAmej74OTxzw1/5NUsKGYLq0bFKYLttaBh+Z6+3y+cXdj4dFOcYo0sST5CyEaiKomfYf+MOB7sPQfcCQ5q3uGGorpX7sn1Hj9wWorZ9q286prLJ0LUnN0UyBJ/kKIBqKuST/mftBumP+HiF4/0WsChCtj7Uv6oS4SN2Z9wnGdyxzb+Sk7tDOQJH8hRANR16Rv3R3OuhNW/RsqykK+dqJXuvLdwA3enXMq6RtdJBx2RffcGi63LWZu1gU8cPXwlB3aGUiSvxCigWAt4JAzWc/7FeS3hbn3exY1DyLcpwozPxX4X2iC8Z+kNbGkiEeuHtCgvMP0awfx2UXl5Cgn19zxYNokfpChnkKIAFPGFzcYGgkR1KTPbQWjfwMf/sqzru3plxnuFupTRTRDMiMRrMVvVwq31obj/RvN+nW74G/PeWY0d+gfdQzJTFr+QogGjFrAEc1kHXILtCuGeQ9A3UnDXUJ9qjB7/dtgFxq31nw77VIWTR0T/nf65n3Pco0jftSkGJKZJH8hRCMTS4pYNHVM5EkSwJ4FE/4Eh7bBl08Y7hJqpSuz179tUveVP609Q1jb9ILTL29SDMlMkr8QwjynXQinX+Ep+nZ4e6OHQ32qiCRZR3NPIOYlFXcsgooVcPbdYLMnfJRSvEmfvxDCXBMegS2fwpyp8P3XGz0crJpmuHsN0d4TiHlJxUVPQH47GPx90+9HJANJ/kIIU/jXxflV82u5e9OLsGE29LskoueHS9ah7gkES8BNXlJxz1rYPBcu+A04Qt+PkOQvhMhYgS3jvx4fy/icTyl675fk9xwFOc0jeh2jZO27qAQbshmXlbQWPgI5LesXqzH7fkQykD5/IUTMAlvGdWRxb+2t5FbvgU8ebPLrRjJW3/SVtMpXwIYPPH39+W1CHiNVV/ECSf5CCBMYJedluh//qpvgWfBl22dNet1ws3PjspLWgj9CXhs468f1m2K+eZyEpNtHCBEzu1K4DGb2/o/rOm5rv4kTb/+Yia7pbDlCVDdeQ3WrFOQ5eOiKM8ztc9+xGLZ8Ahc93GCBmphvHichSf5CiJgZJX6AEzqb//T/Hed+fhM3u57jfm6LaqSMb9F0I81yssxNvm43zP0NNO8IwxovTN/km8dJSrp9hBAxKwrS990638G9y5rxnOtibsz6lLG25UDkM3dDdatEcrM1qrH5q/7tGdd/0cMpuyh9NCT5CyFiNmV8MQ67arT9eE0d5ZXV/KVuEmvdPXjM8TSdOQBElrwnlhTROt9h+Fi4m61RVRCtOQKfPARdhsPASWHjSgeS/IUQMZtYUkSz7Ma9yE63xq4UJ8nmJ86fYsfN37OfIIu6kMnbv8WuNY0uLA674sTJupAt+qhqBS2cBicOwCWPgmp8EUtHkvyFEKY4Uu003O7SmjyHnR26I1Od/48hti3cm/1W0C6dwBZ7ZbUTtKcLSeH5ivZsD9Wij3hs/ndLPYvQD70FOqfH+ryRkOQvhDBFsJa8r35PUUEes91nMcM+gdts7zPRscRwf6MWu9Ot0dpzjMCF1cG4RR/R2PzaKpj5Y2jVFcb9PtyvmFYk+QshTBFqLLx/ldBr7n0Jup0NM++CXcsavU6wFntltTPkZK/A50U0Nv/T38GhrTDxSchpEe5XTCuS/IUQpoh4HYCsHLjuVWjREV67Hiq/a/BwU2fNFgTcGA4bz4bZsPRpGH4H9BzVpGOmMqVDLLlmpdLSUr1sWeNWgRAiTezfCM9eBM3awg/neC4GNK4TFKmCPAcrHxwX2c4HNsM/x0Db3vDDj8CRG230SUsptVxrXRpuP5nkJYSwRmEx3PAWvHwVvDQRfji7vpZOrsNWn/wL8hwoBYerjG8o+wS74QwNK46e1kozw/EALe0OmPRyWiX+aEi3jxDCOt1GwPWveVb/evFy5ixeyb3vrGmQ6E/Wubl0YKdG/feBgnUX+Y8eyqaWh6r+RP6x7Xwx+C9Q0NXUXyeVSPIXQsQsplWuep3vWfTl0LcM/HgSHep2NXi42uliwYb99f33AIEj8UMVWfONHsrGyd8cf2OkfR33OO/k1ysKovkV044kfyFETKKaSRtM7zFwy/vkuquYkf0Q59rWNHi4orK6fsTQ9mmX8vh1gyNeYL6ispo8anjW8RfG2ZfzW+fNzHSfm9K1+M0gff5CiJiYtspV0VDuypnGwzWP8JJjGk+4ruJvdVfhwt6oSyeaImulLY/wcM2f6Kt2McV5B2+5LgBSuxa/GaTlL4SIiZmrXF1/8Wiu13/kXfe5/DzrHd7LfoARjm1Nq5uvNZS9yr/d/01ndYgfOv+7PvGnei1+M0jLXwgRk2Bll5vSsq6vm/9RC+YfK+Gh7Fd4nQdQW8ug/d1QNCT8i2gNO76E+X+A777E0e0clvZ9kK1fVKHSpBa/GST5CyFiMmV8caNx+ZG0rP2HX/oSMngXTDlSAwWjWTrmJi47+m/46llYOwO6DIPTL4du50D700+tDVxz1DNv4NvPYN27sHctNO8Alz0OQ25hnM3GuHPjdgpSkkzyEkLEzCiRh2pZG03kctgVaBrU7clz2D03c09vAStegtVvwJ7Vp14oKxe0G1y1p7Z1HgJDfuApzZzdzNTfMxVEOslLkr8QIuFGTpsfsk6Pv6KCPBZNHXNqw9EKzyLrB7dA1UFQNsgrgLanQdcR0Lx9nKJODTLDVwiRtKK5Gdxo35adPf9ETGIa7aOUaqOUmqeU2uz92jrEvi2VUuVKqb/HckwhROqL5mZwpg/JjJdYh3pOBT7VWvcBPvX+HMzvgc9iPJ4QIg0YlVt22BUOW8O5uzIkM35iTf5XAi96v38RmGi0k1JqKNAB+DjG4wkh0oBRueXp1w5i+vcGRTxzV8Qm1j7/Dlrr3QBa691KqUZ3WpRSNuAx4CZgbIzHE0KkiWCzdCXZJ0bY5K+U+gToaPDQfREe4y5gttZ6pwqzMLJS6g7gDoBu3bpF+PJCCCGiFTb5a60vDPaYUmqvUqqTt9XfCdhnsNvZwHlKqbuA5kC2Uuq41rrR/QGt9TPAM+AZ6hnpLyGEECI6sXb7zAJuBqZ5v74XuIPW+gbf90qpW4BSo8QvhBBNEe0EM+ER6w3facBFSqnNwEXen1FKlSqlno01OCGECMWUctIZSmb4CiFSVrCZwo1mBWeQSGf4SklnIUTKMrOcdKaR5C+ESFkF+Y6ototTJPkLIVJWsF7rJO3NTiqS/IUQKetItTOq7eIUSf5CiJQVrOibFIMLT0o6CyFSiv+4/lZ5Dhx2hdPVcAEYKQYXniR/IUTKCFwBrLLaicOmaJ3voLLKGXKSl0wGa0iSvxAiZUyfu7HB0o/gWfYxPzuLst+OC/q8wIuGbzIYZG4hOenzF0KkjKaO6ze6aFQ7XUyfu9G02FKNJH8hRMpo6g1emQzWmCR/IUTKMFoBLJIbvDIqqDFJ/kKIlGG0Algkq3019aKRzuSGrxAipQRbASzccwAZ7eNHkr8QIiM05aKRzqTbRwghMpAkfyGEyECS/IUQIgNJ8hdCiAwkyV8IITKQJH8hhMhAkvyFECIDyTh/IURKk1LNTSPJXwiRsqRUc9NJt48QImVJqeamk+QvhEhZUqq56ST5CyFSlpRqbjpJ/kKIlCWlmptObvgKIVKWlGpuOkn+QoiUJqWam0a6fYQQIgNJ8hdCiAwkyV8IITKQJH8hhMhAkvyFECIDSfIXQogMJMlfCCEykCR/IYTIQJL8hRAiA0nyF0KIDCTJXwghMpAkfyGEyECS/IUQIgMprbXVMRhSSu0HdsT4Mu2AAyaEYzaJK3LJGBNIXNGSuKITS1zdtdaF4XZK2uRvBqXUMq11qdVxBJK4IpeMMYHEFS2JKzqJiEu6fYQQIgNJ8hdCiAyU7sn/GasDCELiilwyxgQSV7QkrujEPa607vMXQghhLN1b/kIIIQykfPJXSn1PKbVOKeVWSgW9O66UmqCU2qiU2qKUmuq3vadSaqlSarNS6g2lVLYJMbVRSs3zvuY8pVRrg31GK6VW+v2rUUpN9D72glLqW7/HBscaU6Rxefdz+R17lt92089VpHEppQYrpRZ7/9arlVLX+T1m6vkK9l7xezzH+/tv8Z6PHn6P3evdvlEpNT6WOJoQ1y+VUuu95+dTpVR3v8cM/6YJiusWpdR+v+Pf7vfYzd6/+2al1M0Jjutxv5g2KaUq/R6Ly/lSSj2vlNqnlFob5HGllHrCG/NqpdQQv8fMPVda65T+B5wOFAMLgdIg+9iBrUAvIBtYBfT3PvYmMNn7/dPAj02I6VFgqvf7qcCfw+zfBjgE5Ht/fgG4Ng7nKqK4gONBtpt+riKNC+gL9PF+3xnYDRSYfb5CvVf89rkLeNr7/WTgDe/3/b375wA9va9jT2Bco/3eQz/2xRXqb5qguG4B/m7w3DbANu/X1t7vWycqroD9/wt4PgHnaxQwBFgb5PFLgDmAAs4ClsbrXKV8y19r/Y3WemOY3YYDW7TW27TWtcDrwJVKKQWMAd727vciMNGEsK70vlakr3ktMEdrXWXCsUOJNq56cTxXEcWltd6ktd7s/b4C2AeEncjSBIbvlRDxvg2M9Z6fK4HXtdYntdbfAlu8r5eQuLTWC/zeQ0uALiYdO6a4QhgPzNNaH9JaHwbmARMsiut64DWTjh2U1vo/eBp6wVwJvKQ9lgAFSqlOxOFcpXzyj1ARsNPv513ebW2BSq11XcD2WHXQWu8G8H5tH2b/yTR+4/3R+7HvcaVUjgkxRRNXrlJqmVJqia8rividq2jiAkApNRxPa26r32azzlew94rhPt7zcQTP+YnkufGMy99teFqQPkZ/00TGdY337/O2UqprlM+NZ1x4u8d6AvP9NsfrfIUTLG7Tz1VWLE9OFKXUJ0BHg4fu01q/F8lLGGzTIbbHFFMkz/d7nU7AAGCu3+Z7gT14EtwzwK+BhxMYVzetdYVSqhcwXym1BjhqsF/EQ8VMPl8vAzdrrd3ezU0+X0aHMNgW+Hua/n6KQMSvrZS6ESgFzvfb3OhvqrXeavT8OMT1PvCa1vqkUupOPJ+axkT43HjG5TMZeFtr7fLbFq/zFU7C3lspkfy11hfG+BK7gK5+P3cBKvDUzihQSmV5W3C+7THFpJTaq5TqpLXe7U1W+0K81CTgXa210++1d3u/PamU+hdwTyQxmRWXt1sFrfU2pdRCoASYQRPPlVlxKaVaAh8C93s/Evteu8nny0Cw94rRPruUUllAKzwf5SN5bjzjQil1IZ4L6vla65O+7UH+pmYks7Bxaa0P+v34T+DPfs+9IOC5C02IKaK4/EwGfuK/IY7nK5xgcZt+rjKl2+droI/yjFbJxvPHnqU9d1IW4OlzB7gZiOSTRDizvK8VyWs26mv0JkBfP/tEwHBkQDziUkq19nWbKKXaASOB9XE8V5HGlQ28i6c/9K2Ax8w8X4bvlRDxXgvM956fWcBk5RkN1BPoA3wVQyxRxaWUKgH+AVyhtd7nt93wb5rAuDr5/XgF8I33+7nAOG98rYFxNPwEHNe4vLEV47mButhvWzzPVzizgB94R/2cBRzxNm7MP1fxuKOdyH/AVXiuiieBvcBc795ZWVMAAAEDSURBVPbOwGy//S4BNuG5et/nt70Xnv+gW4C3gBwTYmoLfAps9n5t491eCjzrt18PoBywBTx/PrAGTxJ7BWhu0rkKGxdwjvfYq7xfb4vnuYoirhsBJ7DS79/geJwvo/cKnm6kK7zf53p//y3e89HL77n3eZ+3EbjY5Pd6uLg+8f4f8J2fWeH+pgmK6xFgnff4C4B+fs+91XsetwA/TGRc3p8fAqYFPC9u5wtPQ2+39728C8+9mTuBO72PK+BJb8xr8BvBaPa5khm+QgiRgTKl20cIIYQfSf5CCJGBJPkLIUQGkuQvhBAZSJK/EEJkIEn+QgiRgST5CyFEBpLkL4QQGej/Ax2Uo4kMFWRMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x,y)\n",
    "plt.plot(xpl_left,ypl_left)\n",
    "plt.plot(xpl_right,ypl_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 701,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([  62.07662454,   99.98500403, 1442.83746254])"
      ]
     },
     "execution_count": 701,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#recall advertising budget problem on section 12.4:\n",
    "R = np.vstack([[.97,1.86,.41],[1.23,2.18,.53],[.8,1.24,.62],[1.29,.98,.51],[1.1,1.23,.69],[.67,.34,.54],[.87,.26,.62],[1.1,.16,.48],[1.92,.22,.71],[1.29,.12,.62]])\n",
    "m,n = np.shape(R)\n",
    "vdes = 1e3 * np.ones(m)\n",
    "s = npl.lstsq(R,vdes)[0]\n",
    "s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 702,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[315.16818459],\n",
       "       [109.86643348],\n",
       "       [858.96538193]])"
      ]
     },
     "execution_count": 702,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#one can add a total budget constraint given these tools:\n",
    "cls_solve(R,vdes, np.ones((1,n)), np.vstack([1284]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 703,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.05454545],\n",
       "       [ 0.04242424],\n",
       "       [ 0.03030303],\n",
       "       [ 0.01818182],\n",
       "       [ 0.00606061],\n",
       "       [-0.00606061],\n",
       "       [-0.01818182],\n",
       "       [-0.03030303],\n",
       "       [-0.04242424],\n",
       "       [-0.05454545]])"
      ]
     },
     "execution_count": 703,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#minimum norm force sequence\n",
    "A,b = np.eye(10),np.zeros((10,1))\n",
    "C = np.vstack([np.ones((1,10)),np.linspace(9.5,.5,10)])\n",
    "d = np.vstack([0,1])\n",
    "fln = cls_solve(A,b,C,d)\n",
    "fln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.2 Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 704,
   "metadata": {},
   "outputs": [],
   "source": [
    "A,b,C,d = np.random.randn(10,5),np.random.randn(10),np.random.randn(2,5),np.random.randn(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 705,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cls_solve_kkt(A,b,C,d):\n",
    "    m,n = np.shape(A)\n",
    "    p,n = np.shape(C)\n",
    "    G = np.matmul(A.T,A) #Gram\n",
    "    KKT = np.vstack([np.hstack([2*G, C.T]),np.hstack([C,np.zeros((p,p))])])\n",
    "    xzhat = np.vstack(npl.lstsq(KKT,np.hstack([np.matmul(2*A.T,b),d]))[0])\n",
    "    return xzhat[0:n,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 706,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-3.55271368e-15],\n",
       "       [ 3.21964677e-15]])"
      ]
     },
     "execution_count": 706,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = cls_solve_kkt(A,b,C,d)\n",
    "np.matmul(C,x) - np.vstack(d) #residual small"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.3 Solving Constrained Least Squares Problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 707,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.46964193],\n",
       "       [ 0.45536505],\n",
       "       [-0.09061721],\n",
       "       [-0.02080793],\n",
       "       [-0.24730542]])"
      ]
     },
     "execution_count": 707,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.46964193],\n",
       "       [ 0.45536505],\n",
       "       [-0.09061721],\n",
       "       [-0.02080793],\n",
       "       [-0.24730542]])"
      ]
     },
     "execution_count": 707,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def cls_solve(A,b,C,d):\n",
    "    m,n = np.shape(A)\n",
    "    p,n = np.shape(C)\n",
    "    Q,R = npl.qr(np.vstack([A,C]))\n",
    "    Q1 = Q[0:m,:]\n",
    "    Q2 = Q[m:m+p,:]\n",
    "    Qtil, Rtil = npl.qr(Q2.T)\n",
    "    first = np.vstack(np.matmul(np.matmul(2*Qtil.T , Q1.T),b)) #splitting line below for readability\n",
    "    denom = first - np.vstack((2*(npl.lstsq(Rtil.T,d)[0]))) #splitting line below for readability\n",
    "    w = npl.lstsq(Rtil, denom)[0]\n",
    "    return npl.lstsq(R,np.vstack(np.matmul(Q1.T,b)) - np.matmul(Q2.T , w)/2)[0]\n",
    "\n",
    "m,n,p = 10,5,2\n",
    "A,b,C,d = np.random.randn(m,n), np.random.randn(m),np.random.randn(p,n),np.random.randn(p)\n",
    "xKKT = cls_solve_kkt(A,b,C,d)\n",
    "xQR = cls_solve(A,b,C,d)\n",
    "xKKT\n",
    "xQR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 709,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.494448843595619e-16"
      ]
     },
     "execution_count": 709,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npl.norm(xKKT-xQR) #small residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 744,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.6740042981298533e-14"
      ]
     },
     "execution_count": 744,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#sparse constrained least squares\n",
    "def cls_solve_sparse(A,b,C,d):\n",
    "    m,n = np.shape(A)\n",
    "    p,n = np.shape(C)\n",
    "    one = np.hstack([np.zeros((n,n)),A.T, C.T])\n",
    "    two = np.hstack([A, -np.eye(np.shape(A)[0])/2, np.zeros((m,p))]) #2nd object written in J as -I/2, wow\n",
    "    three = np.hstack([C, np.zeros((p,m)), np.zeros((p,p))])\n",
    "    bigA = np.vstack([one,two,three])\n",
    "    xyzhat = npl.lstsq(bigA,np.vstack([np.zeros((n,1)),np.vstack(b),np.vstack(d)]))[0]\n",
    "    xhat = xyzhat[:n]\n",
    "    return xhat\n",
    "    \n",
    "m,n,p = 100,50,10\n",
    "A,b,C,d = np.random.randn(m,n), np.random.randn(m),np.random.randn(p,n),np.random.randn(p)\n",
    "x1 = cls_solve(A,b,C,d)\n",
    "x2 = cls_solve_sparse(A,b,C,d)\n",
    "npl.norm(x1-x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 754,
   "metadata": {},
   "outputs": [],
   "source": [
    "p,n = 50,500\n",
    "C = np.random.randn(p,n)\n",
    "d = np.random.randn(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 766,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  if __name__ == '__main__':\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:10: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n",
      "/Users/vb/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:11: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n",
      "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.01659012],\n",
       "       [-0.00313529],\n",
       "       [ 0.005575  ]])"
      ]
     },
     "execution_count": 766,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.01659012],\n",
       "       [-0.00313529],\n",
       "       [ 0.005575  ]])"
      ]
     },
     "execution_count": 766,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.01659012],\n",
       "       [-0.00313529],\n",
       "       [ 0.005575  ]])"
      ]
     },
     "execution_count": 766,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1 = np.vstack(npl.lstsq(C,d)[0])\n",
    "x2 = cls_solve(np.eye(n), np.zeros((n,1)), np.vstack(C), np.vstack(d))\n",
    "x3 = np.vstack(np.matmul(npl.pinv(C),d))\n",
    "x1[:3]\n",
    "x2[:3]\n",
    "x3[:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 767,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.665155985652124e-15"
      ]
     },
     "execution_count": 767,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "5.691469220174031e-15"
      ]
     },
     "execution_count": 767,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "npl.norm(x1-x2)\n",
    "npl.norm(x2-x3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
